function [ Params , Inversion ] = EstimateRest( options , Params , Targets , data )


%% USEFUL DEFINITIONS

Params.sigma = Params.delta / Params.dhat ;
Params.kappa = ( Params.dhat + sqrt( Params.mu + Params.Delta + Params.dhat^2) ) / Params.sigma ;
Params.tau = Params.delta / Params.sigma^2 * ( sqrt(1 + 2 * Params.rho * Params.sigma^2 / Params.delta^2 ) - 1 ) ; % adimensional
Params.zl0   = Params.tau * ( Params.rho + Params.delta - Params.sigma^2/2) / (1 + Params.tau) ; ...

% assign data
st  = data.st ;
s   = data.s ;
f   = data.f ;
W   = data.W ;
pop = data.pop ;
u   = data.u ;

% recover zeta, v and bpv
run InvertZetaV.m

%% ESTIMATION OF DISTRIBUTION OF ZETA
% Use mean and variance

% Truncate extreme values
[ zetaS, indS ] = sort(zeta) ;
popS = data.pop(indS) ;
cumpop = cumsum( popS ) ;
zetaStrim = zetaS( cumpop < 0.99 )  ;

zetaR = ( zetaStrim - min(zetaStrim) ) / ( max(zetaStrim) - min(zetaStrim) ) ;
mz = mean( zetaR  ) ;
vz = mean( zetaR.^2 ) - mz^2 ;

p2 = mz * ( 1 - mz )^2 / ( 2 - mz ) / vz ;
p1 = mz / ( 1 - mz ) * p2 ;

Params.p1zeta = p1 ;
Params.p2zeta = p2 ;
Params.TailShiftZeta = 0 ;
Params.zetaL = min(zetaStrim) ;
Params.zetaU = max(zetaStrim) ;


%% ESTIMATION OF zhat0

% Quarterly total job finding probability obtained from duration
xPlusy = 1 / Targets.UnempDuration ;
y = ( 1 - Targets.FractionANPE ) * xPlusy ;
Ss = 1 / Targets.QuartersSinceContact - y ;
Ssomega = 1 / Targets.QuartersSinceOffer - y ;
AcceptanceProba = Ssomega / Ss ;

OtC = @( zhat0 ) sum( data.pop .* ( zhat0 * Params.rho / Params.zl0 ./ ( Params.b + v ) ).^zeta ) ...
                  - AcceptanceProba ;

zhat0Grid = linspace(0.001,2,1000)' ;
for i = 1:1000
    otcGrid(i) = OtC(zhat0Grid(i)) ;
end

Params.zhat0 = FindZeroInterp( otcGrid , zhat0Grid ) ;
Params.B0    = Params.rho * Params.zhat0 / Params.zl0 ;


%% HOUSE PRICE ELASTICITY
  
hRes = HousePriceElasticity( data , Params , options  ) ;
Params.eta = hRes.eta ;
    

%% PRODUCTIVITIES

Data.u = data.u ;
Data.s = data.s ;
Data.st = data.st ;
Data.f = data.f ;
Data.W = data.W ;
Data.pop = data.pop ;
Data.houseprice = data.houseprice ;
Data.area = data.area ;

output = inversion( Data , Params , options ) ;
lL = output.lL ;


%% ESTIMATE NU AND AMENITIES

Data0.u    = data.u0 ;
Data0.s    = data.s0 ;
Data0.st   = data.st0 ;
Data0.f    = data.f0 ;
Data0.W    = data.W0 ;
Data0.pop  = data.pop0 ;
Data0.area = data.area ;

Data1.u    = data.u1 ;
Data1.s    = data.s1 ;
Data1.st   = data.st1 ;
Data1.f    = data.f1 ;
Data1.W    = data.W1 ;
Data1.pop  = data.pop1 ;
Data1.area = data.area ;

output0 = inversion( Data0 , Params , options ) ;
output1 = inversion( Data1 , Params , options ) ;

lhs = log( [ data.mig0 , data.mig1 ] ) ;
rhs = log( [ output0.FlowUNetAm , output1.FlowUNetAm ] ) ;
iv  = log( [ data.indemp0 , data.indemp1 ] ) ;

LHS = lhs(:,2) - lhs(:,1) ;
RHS = rhs(:,2) - rhs(:,1) ;
IV  = iv(:,2)  - iv(:,1) ;

drop = (   isnan( LHS ) + isnan( RHS ) + isnan( IV ) ...
         + isinf( LHS ) + isinf( RHS ) + isinf( IV ) ) > 0 ;

% drop too small migration shares for robustness
drop = drop | ( data.pop0 < 0.01 ) ;
     
LHS(drop) = [] ;
RHS(drop) = [] ;
IV(drop) = [] ;

GravityReg = MyTsls( LHS-mean(LHS) , RHS-mean(RHS) , IV-mean(IV)  ) ;
Params.nu  = GravityReg.Coefficients.Estimate(2) ;
Params.eps = 1 / Params.nu ;
Params.P1  = ( 1 + Params.eta + Params.psi ) * Params.eps + Params.omega + Params.psi ;
          
%% AMENITIES AND VARIANCES

Data.HousePrice = output.HousePrice ;
Data.lL = output.lL ;
Data.Ez = output.Ez ;
[ lA , ~ ] = inversionAm( Data , Params  ) ;

% Fit this joint distribution with a lognormal
 
Params.stdlL = sqrt(var(lL)) ;
Params.stdlA = sqrt(var(lA)) ;
Params.corrlLlA = corr(lL,lA) ;
Params.stdx = sqrt(   ( Params.omega + Params.eps )^2 * Params.stdlL^2 ...
                    + Params.psi^2                    * Params.stdlA^2 + ...
                  - 2 * ( Params.omega + Params.eps ) * Params.psi * ...
                        Params.corrlLlA * Params.stdlL * Params.stdlA  ...
                   ) ...
              / Params.P1 ;

%% ASSIGN INVERSION

Inversion.zeta = zeta ;
Inversion.v = v ;
Inversion.bpv = bpv ;
Inversion.lL = output.lL  ;
Inversion.lA = lA  ;
Inversion.HousePrice = output.HousePrice  ;
Inversion.zetaR = zetaR ;
Inversion.Ez = Ez ;

%% P2
Params.P2 = ( Params.omega + Params.eps * ( 1 + Params.eta ) ) / Params.P1 ;



end

